{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# How to debug a model\n",
    "\n",
    "There are various levels on which to debug a model. One of the simplest is to just print out the values that different variables are taking on.\n",
    "\n",
    "Because `PyMC3` uses `Theano` expressions to build the model, and not functions, there is no way to place a `print` statement into a likelihood function. Instead, you can use the `Theano` `Print` operatator. For more information, see:  theano Print operator for this before: http://deeplearning.net/software/theano/tutorial/debug_faq.html#how-do-i-print-an-intermediate-value-in-a-function.\n",
    "\n",
    "Let's build a simple model with just two parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "CompoundStep\n",
      ">Metropolis: [sd]\n",
      ">Metropolis: [mu]\n",
      "Sampling 2 chains: 100%|██████████| 11000/11000 [00:02<00:00, 4276.20draws/s]\n",
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/mkl_fft/_numpy_fft.py:1044: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  output = mkl_fft.rfftn_numpy(a, s, axes)\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import theano.tensor as tt\n",
    "\n",
    "x = np.random.randn(100)\n",
    "\n",
    "with pm.Model() as model:\n",
    "    mu = pm.Normal('mu', mu=0, sigma=1)\n",
    "    sd = pm.Normal('sd', mu=0, sigma=1)\n",
    "\n",
    "    obs = pm.Normal('obs', mu=mu, sigma=sd, observed=x)\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(5000, step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., ..., 0., 0., 0.])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['mu']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hm, looks like something has gone wrong, but what? Let's look at the values getting proposed using the `Print` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Only 3 samples in chain.\n",
      "Sequential sampling (1 chains in 1 job)\n",
      "CompoundStep\n",
      ">Metropolis: [sd]\n",
      ">Metropolis: [mu]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sd __str__ = -1.0187422257064631\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = -0.4591277292273646\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.5731278873643482\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = 1.9855846003705118\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.05734873683282058\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = 1.0169032920056222\n",
      "mu __str__ = 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/twiecki/working/projects/pymc/pymc3/sampling.py:467: UserWarning: The number of samples is too small to check convergence reliably.\n",
      "  warnings.warn(\"The number of samples is too small to check convergence reliably.\")\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    mu = pm.Normal('mu', mu=0, sigma=1)\n",
    "    sd = pm.Normal('sd', mu=0, sigma=1)\n",
    "    \n",
    "    mu_print = tt.printing.Print('mu')(mu)\n",
    "    sd_print = tt.printing.Print('sd')(sd)\n",
    "    \n",
    "    obs = pm.Normal('obs', mu=mu_print, sigma=sd_print, observed=x)\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(3, step, tune=0, chains=1, progressbar=False) # Make sure not to draw too many samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code above, we set the `tune=0, chains=1, progressbar=False` in the `pm.sample`, this is done so that the output is cleaner.\n",
    "\n",
    "Looks like `sd` is always `0` which will cause the logp to go to `-inf`. Of course, we should not have used a prior that has negative mass for `sd` but instead something like a `HalfNormal`.\n",
    "\n",
    "We can also redirect the output to a string buffer and access the proposed values later on (thanks to [Lindley Lentati](https://github.com/LindleyLentati) for providing this example):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO (theano.gof.compilelock): Waiting for existing lock by process '70988' (I am process '70801')\n",
      "INFO (theano.gof.compilelock): To manually release the lock, delete /Users/twiecki/.theano/compiledir_Darwin-18.5.0-x86_64-i386-64bit-i386-3.6.7-64/lock_dir\n",
      "INFO (theano.gof.compilelock): Waiting for existing lock by process '70988' (I am process '70801')\n",
      "INFO (theano.gof.compilelock): To manually release the lock, delete /Users/twiecki/.theano/compiledir_Darwin-18.5.0-x86_64-i386-64bit-i386-3.6.7-64/lock_dir\n",
      "Only 5 samples in chain.\n",
      "Sequential sampling (1 chains in 1 job)\n",
      "CompoundStep\n",
      ">Metropolis: [sd]\n",
      ">Metropolis: [mu]\n",
      "/Users/twiecki/working/projects/pymc/pymc3/sampling.py:467: UserWarning: The number of samples is too small to check convergence reliably.\n",
      "  warnings.warn(\"The number of samples is too small to check convergence reliably.\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAEKCAYAAADdKRa4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvXl4W2eZsH8/XmXHsuMlltPYJU0ap2zpQmibQqG0pZQdZljaGaDMwDAfOxQYhm/4WD/mB5R1LrYBWpZvOmxlK0xLW0ppC0mgLXRv4yxNazex7NiOJceWbFnP74/zHkdR5V3LkfPc16VL0nvOu5yjxM951ldUFcMwDMMwgklFqRdgGIZhGMbsmKA2DMMwjABjgtowDMMwAowJasMwDMMIMCaoDcMwDCPAmKA2DMMwjABjgtowDMMwAowJasMwDMMIMCaoDcMwDCPAVJV6AUb509bWpuvXry/1MgzDMMqKu+6665CqrpnvPBPUxrJZv349d955Z6mXYRiGUVaIyKMLOc9M34ZhGIYRYExQG4ZhGEaAMUFtGIZhGAHGBLVhGIZhBBgT1IZhGIYRYIouqEXkFBHZISJJEXn/HOddLSK7ROR+EblKRKrn6+/OGxCR+7PaW0TkJhHZ7d6bXbuIyH+IyB4RuVdEzsjo8xk39/0i8tqM9pNE5E9urB+JSE3GsdeIyIMi8oCI/HdG+2Xu/N0icllG+6dEpFdExrLWe6KI3CIif3XrepFrrxaR74nIfSLykIh8KKPPe92894vID0Qk5Nrf4a5PRaQt4/y/d2PfKyLbReTUjGOrReQaEXnYzbNttt/JMAzDKCyl0KiHgXcBn5vnvKuBU4CnA3XAmxfQ/7vAxTna/xW4WVU3ATe77wAvBDa511uArwOIyIuBM4DTgLOAD4hIo+vzGeCLbqwR4E2uzybgQ8CzVPWpwHtcewvwUTfOmcBH/QcF4FeuLZsPAz9W1dOBS4CvufZXA7Wq+nTgGcA/i8h6EVnn7slWVX0aUOn6AfwRuBDITgN4BHiuqm4BPgl8M+PYl4HfqOopwKnAQznWaBiGYRSBogtqVR1Q1TuAqXnOu04dwJ+Bzvn6q+pteII8m5cD33Ofvwe8IqP9+26ancBqEVkLPAW4VVVTqnoEuAe4WEQEOB+4JsdY/wR8VVVH/HW69hcAN6nqsDt2E+5hQlV3qurBXJcP+A8GTcCBjPZVIlKF9/AyCcTcsSqgzh2r9/uo6l9VdX+Oe7XdXyuwE3d/3QPJc4Ar3XmTqno4xxoNwzDKit/c3080lij1MhZN4H3UzuT9euA3yxgm4gtE997u2tcBvRnn9bm2e4AXiki9Mxc/D+gCWoHDqprKOh+gG+gWkT+KyE4R8TX72eaYi48BrxORPuA64J2u/RrgCHAQeAz4nHsAeBzPwvCYOzaqqjfOM0cmbwKud583AIPAd5zp/dsisiq7g4i8RUTuFJE7BwcHFzGVYRhG8ZmYnOatV9/F/9uxoBojgSLwghrP7Hubqt5egLElR5s6IXcdsB34AbADSM12vnuvwjOhnwdcCnxbRFbP02c2LgW+q6qdwIuA/yciFXhm8mngBOAk4H0issGZ0l/u2k7A07pfN88cAIjI8/AE9QczruMM4OvO9H6Eo66Coxeg+k1V3aqqW9esmbcCnmEYRknpjyVQhYOjplHnRETeLiJ3u9cJi+j3UWANcPkylxB1Jm3cu2+W7sPTlH06OWoy/pSqnqaqz8cTtruBQ3jm8ars891Yv1TVKVV9BNiFJ7hnnWMO3gT82K1jBxAC2oC/w/MdTznT+h+BrXg+6EdUdVBVp4CfAefMd1NEZAvwbeDlqjqUcR19qvon9/0aPMFtGIZRtvgm74G4CeqcqOpXndA7TVXnE1IAiMib8fy7l6pqeplLuBbwo60vA36Z0f4GF/19Np7J+KCIVIpIq1vHFmALcKPzl98CvCrHWL/AM5HjzOXdwD7gBuAiEWl2mu9Frm0uHgMucGM9GU9QD7r28916VwFnAw+79rOdqV5c3zkDwETkRDyB/npV7fHbVbUf6BWRza7pAuDBedZrGIYRaHxBbT7qBSAiHc73ejnwYRHp8yOqReS6DI37G0AE2OE08Y8soL9vpt7s2t/kxvo08HwR2Q08330Hz7y9D9gDfAt4m2uvBm4XkQfxoqFfl+GX/iBwuYjswfNZX+nabwCGXJ9bgA+o6pCqDuNFVd/hXp9wbYjIZ9211Lv1fsyN9T7gn0TkHjzT+xvdQ8JXgQbgfjfWd1T1Xqf9XgP8BbgP73f9ppvjXW6OTuBeEfm2m+Mjbv1fc/c3c1eNdwJXi8i9eJHv/57zxzQMwygTBmJJAKLuvZwQ7++/YSydrVu3qu2eZRhGkPnkrx/kyj88AsDDn7yYUHVliVcEInKXqm6d77xyCCYzDMMwjGWRafIuN/O3CWrDMAxjxTMQS1JZ4SXhlJv52wS1YRiGseLpjyXojoRnPpcTJqgNwzCMFY2qEo0lOLWzCYABE9SGYRiGERxiEymSqTQntzcQqq4wH7VhGIZhBImoK3ISaQwRaQyZj9owDMMwgkT/aIagDofMR20YhmEYQcI3dXc0hog0hcxHbRiGYRhBYiDumbrbG2uJhGuJxpKUU7EvE9SGYRjGiqZ/NEFTXTWh6koijSEmpqaJJVLzdwwIJqgNwzCMFU00lqCjMQRApMl7LyfztwlqwzAMY0UTjSdpb6wFIBL23ssp8tsEtWEYhrGiGYgliPgatXsvp1xqE9SGYRjGimU6rQzEk0R8jdoJ6nJK0TJBbRiGYaxYho4kmU7rjI+6rqaSxlCV+agNwzAMIwgMxPzUrNBMW7lVJzNBbRiGYaxYMquS+UQaQzNlRcsBE9SGYRjGiuVone/ambb2xlqioyaoDcMwDKPkRGNJRGBNw1FB3dEYYiCeJJ0uj+pkJqgNwzCMFctALEFbQy1VlUfFXaQxRCqtDI9PlnBlC8cEtWEYhrFi6Y8ljjF7w1EzeH+ZmL9NUBuGYRgrlmgsOZOa5eMHlg2USUCZCWrDMAxjxTIQSxyTmgWZ1cnKI0XLBLVhGIaxIplMpRk6MkkkfKygXjNT79s0asMwDMMoGQM5UrMAqisraGuoMUFtGIZhGKXEN237W1tmUk7VyUxQG4ZhGCsSv553tukbfEFtGrVhGIZhlAx/h6xs07ffZoLaMAzDMEpINJakulJoWVXzhGORxhCHxiaZmk6XYGWLwwS1YRiGsSIZiCVoD4cQkScc81O0BuPB91MHSlCLyFUiMiAi989y/OUicq+I3C0id4rIszOOTbv2u0Xk2oz280XkLyJyv4h8T0SqXHuTiPxKRO4RkQdE5B8y+lwmIrvd67KM9te6+R8Qkc9mtF8uIg+6YzeLyJOWsa65rvEz7vz7ReS1Ge0iIp8SkR4ReUhE3pXR/h8isseNeYZrf17Gmu4WkYSIvMIduz2j/YCI/GKxv6NhGEYQiMafWJXMx28vB/N3VakXkMV3ga8A35/l+M3AtaqqIrIF+DFwijs2oaqnZZ4sIhXA94ALVLVHRD4BXAZcCbwdeFBVXyoia4BdInI10AB8FNgKKHCXE7AVwBXAM1R10AnXC1T1ZuCvwFZVHReRtwKfBXxButh15bxGEXkxcAZwGlAL3Coi16tqDHgj0AWcoqppEWl3U70Q2OReZwFfB85S1VvcOIhIC7AHuBFAVc/NWOdPgV/O8lsYhmEEmv7RBN2RcM5j7WG/6EnwBXWgNGpVvQ0YnuP4mKr6252swhOkc9EKJFW1x32/CfhbfzggLJ5NpMHNmwJeANykqsOqOuL6XAxsAHpUddD1/60/lqreoqrjrn0n0LnUdc1xjU8BblXVlKoeAe5x6wJ4K/AJVU27MQZc+8uB76vHTmC1iKzNWsurgOsz1g+AiISB8wHTqA3DKEsGYslj9qHOpKOpfKqTBUpQLwQReaWIPAz8D/CPGYdCzlS80zfjAoeAahHZ6r6/Ck/zBE9zfzJwALgPeLcTdOuA3oxx+1zbHjzNdr0zU78iY6xM3gRcv4x1zXaN9wAvFJF6EWkDnpfRZyPwWjfP9SKyybXPdi2ZXAL8IMd1vBK42WnshmGsMN76X3dx5R8eKfUyCsaRZIp4MjWroG6pr6GqQkyjLgSq+nNVPQVPUH4y49CJqroV+DvgSyKy0WmmlwBfFJE/A3E8rRk8zflu4AQ8M/BXRKQReGLUAajTrt8K/Ai4HdifMRYAIvI6PJP5FctYV85rVNUbgeuA7XiCdUdGn1og4eb5FnCVv6Rc15Kx3rXA04Ebcpx3KbkFuN/3Le7B4M7BwcHZTjMMI4BMp5XfPhTl1p6V+383OkdqFkBFhdAerjWNupA4M/lGp12iqgfc+z7g98Dp7vsOVT1XVc8EbgN2uyH+AfiZMwvvAR7B83f3caym3ImndaOqv1LVs1R1G7ArYyxE5ELg34CXqerML7+Edc11jZ9S1dNU9fl4Qtjv0wf81H3+ObAloz3ntTheA/xcVacy5xWRVuBMPI0+J6r6TVXdqqpb16xZM9tphmEEkGgswdS00jc8Pv/JZYovgLN3zsqkvUyKnpSVoBaRk51PGRfBXAMMiUiziNS69jbgWcCD7nu7e68FPgh8ww33GHCBOxYBNgP78LTLi9yYzcBFri1zrGbgbcC33ffTgf/EE9K+f5ilrGuOa6x0AhQXZLYFFwCG50c+331+LuD7vq8F3uCiv88GRlX1YMYtnU1rfjXwa1UN/r9gwzAWTd/IhPd+eIJ0er5Qn/LEr/OdvXNWJh1lIqgDFfUtIj8AzgPaRKQPL/q6GkBVv4EXcPUGEZkCJoDXuujoJwP/KSJpvIePT6vqg27YD4jIS1z711X1d679k8B3ReQ+PO30g6p6yK3jk8Ad7rxPqKof4PZlETk1o90XiFfgBaT9xMnYx1T1ZXg+8MWua7ZrrAZud+PHgNepqm/6/jRwtYi8FxgD3uzarwNehOdfH8ezIvj3ej2etn1rjp/iEjemYRgrkF6nSU+m0gyOzR5wVc7MZ/r2j23fe6hYS1oygRLUqnrpPMc/A3wmR/t2PF9rrj4fAD6Qo/0Anracq89VHPXzzrs+Vb1wlvalrGu2a0zgRX7nGusw8OIc7YqXhparz36eGFjmHzsvV7thGCuD3pGjJu/e4fEVKaj7R5PU11TSUDu7mGtvDBFLpJiYnKauprKIq1scZWX6NgzDMJZP7/AEFS7UNFNorySi8QQdjbmrkvn4/uugm79NUBuGYRxn9I6M89QTmrzPwxMlXk1hGIglaJ/D7A1Hy4iaoDYMwzACRd/wOJvaG1gTrp3xV680onMUO/GZKSMa8HrfJqgNwzCOIyZTaQ7GEnS21NPVXLciTd+qSn8sMWdqFhyNCI+OmkZtGIZhBISDoxOoQldzHV0t9TOpWiuJ0YkpJlPpOVOzABpDVdRVV5rp2zAMwwgOvk+6q6WeruZ6Do4mSJXBnsyLwS92MldqFoCIEGmsNdO3YRiGERx8U3dXSz1dLXVMp5WDATf9Lpb+mRzq+dPO2htDZvo2DMMwgkPv8DhVFUJHY4jO5vqZtpWEb8qez0ftnxONm6A2DMMwAkLvyAQnrK6jskLo8gX1CgsoG3CCek14btM3eObxaCzB0d2Fg4cJasMwjOOI3uFxulrqAFi7OkSFrLxc6mgsyer6akLV81cbizSGSEyliSVS855bKkxQG4ZhHEf0jUzMaNLVlRWsbaqjb4Vp1P2xBJHwwsqitpdB0RMT1IZhGMcJE5PTHBpL0tVSP9PW1VJH7wpL0RqIJYg0LUxQl0MZURPUhmEYxwm+5tzZXDfT1tVcvwKDyZJEFuCfhozqZLHgpmiZoDYMwzhO6J0R1Ec16s7megbiSRJT06VaVl6ZTuuitu5sD5tGbRiGYQSEo8VOMjRq93mlVCgbGksyndYFm77raippDFWZoDYMwzBKT+/wOKHqCtY0HDUL+/7qlZKiNVOVbIGmb4COppAJasMwDKP09I6M09lcf8wezX4EeN8K8VNHF1GVzCfSGDIftWEYhlF6vNSsumPa2sO11FRVrBjT92LKh/q0h02jNgzDMAKAV+yk/pi2igqhc/XK2e5yIJagQqCtoWbBfTqaahmIJ0mng1mdzAS1YRjGccDoxBSxRGrG1J1JZ0v9iqlOFo0laWuopapy4eIt0hhiOq0MHZks4MqWjglqwzCM4wA/V7ozy/Ttt60Ujbo/lliU2RuCn6JlgtowDOM4oC9je8tsuprrOTw+RTwxVexl5Z3oEgT10aInJqgNwzCMEjGTQ53D9O3nUq8E8/dAPDkjeBdKR5OvUQcz8tsEtWEYxnFA38g44VAVTfXVTzg2k6JV5ubvZGqa4SOTi9ao2xpqETGN2jAMwyghvRm7ZmVztOhJeWvUA36xk0Vq1NWVFbSuqjVBbRiGYZSOzH2os2mur2ZVTWXZb84xEF98DrVPR5MJasMwDKNEqCp9IxPHbMaRiYjQ2Vxf9qbvmfKhSxDUkXBwq5OZoDYMw1jhHBqbZGJq+glVyTLpaqkr+2CypZQP9WlvDM1o5EHDBLVhGMYKp3eO1CyfzuZ6ekfGUQ1mda6F0B9LUFNZQXOOgLn5iDTWcmhskslUugArWx4mqA3DMFY4vu95LkHd1VLP+KQXNV2uDMSStDfWHrPpyELpcFr44FjwzN8mqA3DMFY4/oYbuaqS+fhm8XLenGMpxU58/H5BDCgrqKAWkYtFZJeI7BGRf81x/HIReVBE7hWRm0XkSRnHfiMih0Xk11l9rnZj3i8iV4lItWv/gIjc7V73i8i0iLSISJeI3CIiD4nIAyLy7oyxThWRHSJyn4j8SkQas+Y6UUTGROT9GW2rReQaEXnYjbnNtX9MRB7PWMOLXHu1iHzPzfGQiHwoY6z3ujXdLyI/EJGQaz9JRP4kIrtF5EciUrOA+3WiiNzo5nhQRNa7dhGRT4lIjzv2roz2/3C/zb0icsZ8YxmGUZ70jYzT1lBDfU3VrOeshH2pvfKhi0vN8ml3/QaOJ0EtIpXAV4EXAk8BLhWRp2Sd9ldgq6puAa4BPptx7Arg9TmGvho4BXg6UAe8GUBVr1DV01T1NOBDwK2qOgykgPep6pOBs4G3Z6zj28C/qurTgZ8DH8ia64vA9VltXwZ+o6qnAKcCD2We769BVa9zba8Gat0czwD+WUTWi8g64F3u+p8GVAKXuD6fcWNtAkaANy3gfn0fuMJd55nAgGt/I9AFnOKO/dC1vxDY5F5vAb6+gLEMwyhDeodnj/j28bXtcg4oG4gll61R94+WoaAWkW6nvd3vvm8RkQ8vYOwzgT2quk9VJ/EExMszT1DVW1TVf3zbCXRmHLsZiGcPqqrXqQP4c2afDC4FfuDOP6iqf3Gf43iCdZ07bzNwm/t8E/C3Gdf9CmAf8EBGWyPwHOBKN96kqh6e5z4osEpEqvAeLCaBmDtWBdS5Y/XAAfGcK+fjCWKA7wGvcPPlvF/uwaNKVW9y541lnPdW4BOqmnbHfKH7cuD77lbuBFaLyNp5xjKMFcmu/jg90Sf8uVkx9I6Mz2n2BgiHqlldX122GvVYMsVYMrVkQd1SX0N1pRCNl6eP+lt4GuoUgKrey1HNby7WAb0Z3/s4KiBz8SaeqL3OijN5vx74TVZ7PXAx8NMcfdYDpwN/ck33Ay9zn1+Np3kiIquADwIfzxpiAzAIfEdE/ioi33bn+rzDmZGvEpFm13YNcAQ4CDwGfE5Vh1X1ceBzru0gMKqqNwKtwGFVTbn+s923zPvVDRwWkZ+5dV3hLBoAG4HXisidInK9iGxy7bP9PnONlXkv3+LGvHNwcDDH8gyjfPiXa+7hgz+9t9TLKAjTaeXA4Yk5A8l8uprry7boycBMatbSTN8VFUJ7OFS2Pup6Vf1zVlsq55nHkivsLmfcv4i8DtiKZ+5eKF8DblPV27PaXwr80Zm9M+dowBPe71FVX6P9RzxT+F1AGE/bBU9Af1FVx7LGrgLOAL6uqqfjCWDf9/51PKF4Gp7g/bxrPxOYBk4ATgLeJyIbnCB/uWs7AU/rfh0LuG857lcVcC7wfuCZeA8Ub3THaoGEqm7Fe+i6yh9mlnnmGuvoiarfVNWtqrp1zZo1OYYyjPJgOq3sisbZ1R8nnS7f1KTZ6I8lmJrWWcuHZtLVUle2wWT9y8ih9mlvDGZ1soUI6kMishEnLETkVXiCaD76cBqqoxM4kH2SiFwI/BvwMlVdkM1BRD4KrAEuz3H4EpzZO+P8ajwhfbWq/sxvV9WHVfUiVX2G67PXHToL+KyI7AfeA/xvEXmHu6Y+VfU18mvwBDeqGlXVaWdi/haegAb4Ozyf9pQzO/8RT8heCDyiqoOqOgX8DDgHOIRnhvajPo65b7Pcrz7gr87NkAJ+4a/LHfOtCz8HtmS05/p95hrLMFYcvcPjJKbSjE9O8/jh8hRSc9E3k5o1t+kbPI368ZGJsnxgGVhGVTKfjsZgVidbiKB+O/CfwCki8jie4HrrAvrdAWxyEcw1eAL02swTROR0N/bLMnyncyIibwZeAFzq+10zjjUBzwV+mdEmeD7lh1T1C1nnt7v3CuDDwDcAVPVcVV2vquuBLwH/rqpfUdV+oFdENrshLgAedGOszRj6lXhmdfBM2+e7KOtVeAFtD7v2s0Wk3q3xArdGBW4BXuX6X+Zfzxz36w6gWUR81fZ8f114gvZ89/m5QI/7fC3wBreus/FM7wfnGcswVhy7MnzTu/pXnp/a32hjIRp1Z0s9k9NpBgLop52P5VQl84k0lqnp22lWF+JpsKeo6rNVdf8C+qWAdwA34AVw/VhVHxCRT4iI7xe+AmgAfiJeStOMIBeR24GfABeISJ+IvMAd+gYQAXa4Ph/JmPaVwI2qeiSj7Vl4vuzzJSt1Ci8SvQdPcB4AvjPfdQHvBK4WkXvxzNz/7to/K14K1r3A84D3uvavumu8H08IfkdV73Va+TXAX4D78H6Lb7o+HwQuF5E9eD7rK+e6X6o6jWeqvllE7sMza3/L9fk08Leu/f/DRckD1+EFy+1x575tAWMZxoqjJ0M49wysQEE9PI4InLB6IRq1i/wuw4Cy/liCVTWVNNTOnoI2H+2NtcQTKcYnF+LdLR7zXlGWIJyp+KKqn5ivr0tRui6r7SMZny+co++5s7TPumZV/S7w3ay2P5DbH4uqfhkv3WpWVPVjWd/vxjNdZ5+XK5UM5+d+9SzHPgp8NEf7Po6azjPb57pfN3HUrJ3Zfhh4cY52xbOWLHgsw1iJ9AyM0dVSRzp9rNBeKfSOjNPRGKKman4Dqp/C1Ts8zjPXtxR6aXllIJYk0rR0bRq8jTnA29zjpLalC/x8s5CVZGqnIeAlHJs7bBiGUbb09MfZHAm7oLLs+NHyp2949n2osynnXOpoLDEjaJdKR9PR6mQnta2a5+ziMa+gVtXPZ34Xkc+R5Ws2DMMoRyZTafYOjnH+k9tJp5U/7hkiNZ2mqnLlVFfuHRln28bWBZ0bqq6kPVxblqbvaDzBM05snv/EOfBTu4Lmp17Kv8Z6vJQdwzCMsmb/0BFSaWVzJEx3JMzkdJr9Q+UnpGYjmZqmP5ZYsEYNXinRcsulVlWiy6hK5tMe0HrfC/FR38fRPN5KvKCyef3ThmEYQceP8u6OhEm77R13R+Oc3N5QymXljYOHE6jOvWtWNl3Nddyxf6SAq8o/h8enmEylly2ow7VV1NdUBi5FayE+6pdkfE4B0YyqWYZhGGXL7micygphwxrPHynipWu98Olr5+lZHszsQz1P+dBMulrqufaeA0xNp6kuExdANL781CzwgqWDmKI1q6AWET/kLzsMslFEyK78ZRiGUW7sisZZ31pPqNqrkru+ddWKqvntB4UtTqOuJ62eNn5i68L7lZLoTLGTpZUPzaQ9XDtTPCUozKVR34Vn8p6t1KT5qQ3DKGt6omOc0hGe+b6pvWFFFT3pHRmnulIWpWl2ZuRSl42gHs2PRu2PcXfvfHstFZe5cpJPKuZCDMMwikliapr9Q0d42aknzLRt7ghz88MDJFPT1FY9YS+asqN3eJwTVtdRWZGzlEROZvalLqOAMt9U3Z4HjbqjKUT0gQSqOlM3pNQsKKPbbSCxCS+PGgBVvW32HoZhGMFmz8AYqp5w9ul2+dT7Bo/w5LWNJVxdfugdWXgOtc/aphCVFVJWKVrReILm+uq8PFy1h2tJptLEJlI01VfnYXXLZyH7Ub8Zb8/mG/B2lboB+Fhhl2UYhlFYfF90d+SooPaF9krxU/cNjy9oM45MqiorWNsUKquiJ/2jy0/N8vHH8QPUgsBCQvrejbfd4aOq+jy8/ZxtA2LDMMqaXdE4NZUVrM/ww65vXUVVhawIP/X4ZIqhI5MzZUEXQ1dzPX1lpFEPxBN5F9T9o+UlqBOqmgAQkVpVfRjYPE8fwzCMQNPTH2fDmlXHVCGrqapgw5qVEfnt7yu9mIhvn66Wupldt8qBaCyRl4hv8La69McMCgvxUfeJyGq87RJvEpERcuwrbRiGUU70RMfYuv6JJSe7I2Hu6QtW1O9S8IPBFpND7dPVXM9gPElianomdS2oTKeVwXj+TN9+QFqQtvpcyDaXr1TVw24Xqf+Dt+XiKwq9MMMwjEIRT0zx+OGJY/zTPpsjYXqHJwK31eFi8QX1Ukzfnc6vXQ7m70NjSdKan9Qs8OqdN9VVl5fpW0S+LCLnAKjqrap6rapOFn5phmEYhWH3gLdL1uYcgrrbBZTtLvOdtHpHJqirrqStoWbRfbtmtrsMvvnbN1HnS1CDZ/4Okul7IT7qvwAfFpE9InKFiDxhL2bDMIxyoqf/iRHfPn7brjL3U/cOj9PZXLekXOCZXOoy0KjzWZXMp72xlmiZmb6/p6ovAs4EeoDPiMjugq/MMAyjQOyKxqmrrpypwpXJiS311FZVzAjzcqVvZGJJgWQAaxpqqamqmAlICzL9BdCoI40hBspMo/Y5GTgFWA88XJDVGIZhFIGeaJzuSAMVOSp2VVYImyIN5a9Rj4wvKZAMoKJC6GyuK4vqZAOxBBUCbQ3506gjjbUMxJNMp3X+k4vAQnzUvgb9CeB+4Bmq+tKCr8wwDKPDpLbIAAAgAElEQVRA7Oofy2n29umOhMvaRz06PkU8kVqyRg2en7o8TN8J1oRrF1UmdT46GkNMp5WhI8Ewfy9Eo34E2KaqF6vqd1S1/PMWDMM4bhk+MsmhseQxpUOz2RwJ0x9LMDo+VcSV5Q9fwOYy7S8UT6MOvuk7GstfapZPuxsvKLtoLcRH/Q1VPVSMxRiGYRQav5jJpnk0aoCegfI0fy8nNcunq6We0YkpYolgP6x4xU7yK6iDVp2sPHYFNwzDyBO+oM6VmuXjp2iVaylRX6Nerukbgr+LVj6rkvl0BKzetwlqwzCOK3b1x2kMVc35x/2EphANtVVlW0q0d3iCxlAVTXVL3/3J38wjyObvZGqakfEpIuH8atRtDTWIHE39KjXzlhAVkRNztavqY/lfjmEYRmHpicbZ3BGeM79YROiONJStoO4bGV+WNg1HNeogVycbmMmhzq+grqqsoK2hNjApWgup9f0/gAKCtx/1ScAu4KkFXJdhGEbeUVV6omO8ZMvaec/d3BHmN/f3o6pLKhpSSnpHJjh5TcOyxlhdX01DbVWgc6lnqpI15VdQg5ei1R8QQb2QYLKnq+oW974Jr/DJHwq/NMMwjPwyEE8yOjE1Z2qWz6b2MCPjUxwaK6+KyarqNOqlR3yDZ1UIei51IaqS+XhlRINh+l60j1pV/4K3P7VhGEZZsWuO0qHZ+Olb5Wb+HhxLkphKLyvi26cz4LnUMxp1nn3U4KVolY3pW0Quz/haAZwBDBZsRYZhGAXCF7rdkfnNwjM1v/vjPOvktoKuK5/4wV/L1aj9Mf6451Bgzf/RWIKaqgpW1y89aG42IuEQQ0cmmUylqakqbdz1QmYPZ7xq8XzWLy/kogzDMArBrv44bQ21tC6g3GRbQw0tq2rYXWa51H7wV1ceNOqu5nompqYZOhJM87+fmlWIh4jIzL7Updeq59WoVfXjxViIYRhGoekZGGNzx8KCrPzI73LLpfaDv/Jh+vYjx/tGJvJaSztfRGPJgpi94WiAWjSWzMu9XA4LqfW9VUR+LiJ/EZF7/VehFiQiF4vILret5r/mOP5FEbnbvXpE5HDW8UYReVxEvpKj77Uicn/G99NEZKcb604ROdO1/33GtW4XkVMz+lwlIgOZ47j2H2Wsa7+I3J1xbIuI7BCRB0TkPhEJufbXujkeEJHPLuQaReQ3InJYRH6dNf/tGX0OiMgvso4/U0SmReRV890vN8c9bl3fEJHK7HtpGOVGOq3sjsbZ1D6/f9qnOxKmJzqGajA2Z1gIvcPjtDXUUlez/P+2R3Opg+mnjsbzX5XMx38ACIKfeiHpWVcDHwDuA9KFXIwTCF8Fng/0AXeIyLWq+qB/jqq+N+P8dwKnZw3zSeDWHGP/DZBdZf+zwMdV9XoReZH7fh5effPnquqIiLwQ+CZwluvzXeArwPczB1LV12bM9Xlg1H2uAv4LeL2q3iMircCUe78Cb5OTQRH5nohcoKo3z3ONVwD1wD9nzX9uRp+fAr/M+F4JfAa4Ifu+kPt+vUZVY+LZk64BXg38MEdfwygbHj88wfjk9Jw1vrPpjoQZS6Y4MJpg3erl+3yLQW8eIr59fE0yqAFl0dEE53W3F2Rs3/QdhBSthfioB1X1WlV9RFUf9V8FWs+ZwB5V3aeqk3jCYS5/+KXAD/wvIvIMIALcmHmSiDQAlwP/N6u/Ao3ucxNwAEBVt6vqiGvfCXTOdFC9DRiebUFOuL0mY10XAfeq6j2u/5CqTgMbgB5V9QPzfgv87XzXqKo3A7Pa4kQkDJwPZGrU7wR+CgxknZvzfqlqzH2sAmrw7pNhlDWLifj2mYn8LiPzd+/wRN5MtQ21VTTXVweyOtlYMsWRyemCpGYBtKyqobpSApGitRBB/VER+baIXCoif+O/CrSedUBvxvc+1/YERORJeMVXfue+VwCfx9P+s/mkO5b9WPge4AoR6QU+B3woR983Adcv/BI4F4iq6m73vRtQEbnBuQ/+xbXvAU4RkfVO634F0DXXNS6QVwI3+8JWRNa5tm9kjT3X/UJEbsAT7HE8rdo4Tvje9v28/yf3lHoZeWfXIiK+fbrbyytFazqtHDg8seR9qHPR1VIfyOpkM6lZBTJ9iwjt4WCkaC1EUP8DcBpwMfBS93pJgdaTK3RvNm3uEuAap50CvA24TlUzBT0ichpwsqr+PMcYbwXeq6pdwHuBK7P6Pg9PUH9w4ZdwrAaMp5U+G/h79/5KZ+IecfP/CLgd2A+k5rnGpcz/JeCDOcbIeb98VPUFwFq8SP/zs4+LyFucX//OwUHL1ltJ/PefHuOnf+ljJKCRvktldzTOutV1hEMLT+Vpqq+mozE0I+SDzsHRCVJpXXb50Ey6musD6aOOup2t2gukUYNn/g7CxhwL8VGfqqpPL/hKPPo4VqvsxJmjc3AJ8PaM79uAc0XkbUADUCMiY8CjwDNEZD/e9baLyO9V9TzgMuDdrv9PgG/7g4nIFvf9hao6tJDFO834b4BnZF3Trf5WoSJyHV4u+s2q+ivgV679LUC2MM2+xvnmb8VzH7wyo3kr8EPPIk8b8CIRSTHL/VLVmQA+VU2IyLV47oebMudS1W/i+e7ZunWrmcZXCIfGkjNC6U+PDHHx0+YvtVku7IqOsWkR2rTPpjKq+T2TQ53HKOXOljpufLCf6bRSWRGcXGpfgHYUSKMGT1sPwm+/EI16p4g8peAr8bgD2CQiJ4lIDZ6gujb7JBHZDDQDO/w2Vf17VT1RVdcD7we+r6r/qqpfV9UTXPuz8fzC57luB4Dnus/nA7vd+CcCP8MLAOtZxPovBB5W1b6MthuALSJS7wT5c4EH3Tzt7r0ZT8PNfFB4wjUugFcDv1bVmUdAVT1JVde7678GeJuq/mK2+yUiDSKy1q2hCngR8PAi1mCUMTv3HX0m3bF3Qc+nZUFqOs3egbE5t7acjc2RMLujY0yng/88OpNDnadgMvCE/tS0BiKfOBPfd9xeYEE9EAAf9UI06mcDl4nII0ASzzytqrol34tR1ZSIvANPuFUCV6nqAyLyCeBOVfWF9qXAD3X5ORP/BHzZCaQE8BbX/hGgFfia00RTqroVQER+gBcZ3iYifcBHVdU3mV/CsWZnXOT4F/AeQhTP3Pw/7vCXM1K/PpH1UJDzGkXkduAUoMHN/yZV9aO5LwE+vbRbMcMq4FoRqcX7DX5Hln/bWLls3ztEQ20Vp3Y1sX0FCer9Q+NMTqcXFUjm090RJplK0zs8zvq2VQVYXf7oHZmgQuCEPEao+2b03uEJ1jYFJ/I9GkvQUFtFQ+1CxNjSiDSGiCdTHEmmWFXAeeZjITNfXPBVZKCq1wHXZbV9JOv7x+YZ47t4aVTZ7fuBp2V8/wPHmqn99jcDb55l7EvnmPeNs7T/F16K1mLG+tgs7efmanfHzpvt2Dzr+y7ufqlqFKvlftyyY+8QZ53UwjNPauHT1z/MQDxBe4EKShQT33y5mNQsH18L3xWNB15Q9w2Ps7apjurK/JW89APTeofHOfOklryNu1z8qmSFxB8/GkuwYZm7kS2Hheye9WiuVzEWZxhG8Tg4OsEjh46wbWMr52xsBWDnvlkzEcuKnmgcETi5fWk+aiiPFK3ekXHW5THiG45q50HLpY7GkgWL+Pbx/d+lTtEqbaVxwzACg++T3raxlaee0EQ4VMWOvYdKvKr80BON86SWekLVi6/WVV9TRVdLXVlEfvcOT+Q1kAwgVF1JpLE2cLnUnkZdWEHt+79L7Z8vndHdMIxAsX3vEKvrq3lyRyMVFcJZJ7WuGD/1rv74kvzTPpsj4UBE/85FMjVNNJ7IayCZT1fAtrtUVQZiyYKmZsGxpu9SYhq1YRioKjv2DrFtQysVLgXnnI2tPDo0zuOHg6VJLZbE1DT7h8aX5J/26Y6E2Td4hMlUQasoL4sDhxOo5jc1y6erpZ7HR4Lz72BkfIrJ6XRBU7PAq8xWX1NJ/6iZvg3DKDG9wxM8fniCbc43Dcx8Lvc0rX2DR5hO6/I06o4wqbSyf+hIHleWX/yiJPksduLT1VzHwdEJpqaD8aBS6KpkPiJCR2Oo5EVPTFAbhsF254s+J0NQb46EaVlVM3OsXPH3k16uRg0EesvL3gLkUPt0ttSTVjgQEOvKUUFd+K032xtrS15G1AS1YRjs2DfEmnAtGzNSUCoqhG0bWtmxd6istnnMZld/nKoKYX3r0lOrNqxZRWWFBNpP3Ts8QXWlFCSdrnMmRStogrrwqYORxpBFfRuGUVpUle3OP+0K/MywbWMrB0cTPDoUnECixdITjbNhzSpqqpb+5662qpL1rfWB16jXra4rSJnProBtd+kLzjXhwmvUkcYQ/bFESR9WTVAbxnHO3sExBuPJY8zePr6fupyjv3dFlxfx7bO5I8zugewt7YND3/B4QfzTAGubQlRWSGA254jGErSsqqG2avHpdosl0hhiMpVmdGKq4HPNhglqwzjO8YPFztnY9oRjG9pWEWmsLVs/9ZFkit7hiSXV+M6mOxJm/9ARElOL2cyuePSO5G8f6myqKis4YXWI3oBEfkdjSdqLoE1DZopW6czfJqgN4zhn+94h1q2uyxmEJCKcs7GNnfvK00+9x2nA3csIJPPZHAmjenTMIHEkmWL4yGRBAsl8upqDsy91NJago6k4pW19P3h/CQPKTFAbxnFMOq3s2DfEto1P9E/7bNvQyqGxyUCbfWfDryaWD9P3pgBHfveN5H97y2y8famDolEniBSpBv3RMqImqA3DKAEP9cc4PD6V0z/tM+On3lN+5u+e/ji1VRWcmAff7frWemoqKwIZ+V3IHGqfrpY6Do0lmZgsrek/NZ3m0FiyKKlZcDRgrZQpWiaoDeM4JrO+92x0tdTT1VLHjn3lF1C2KxpnU6QhL5HQVZUVbGxvCKagdibpzjxvyJGJ7/8utfn70NgkaYVIkUzfoepKVtdXm4/aMIzSsGPvECe1rZp3n+FtG1rZuW+Y6XR5+al78hTx7bM50kBPNHgugN7hCeqqK2ldVVOwOXz/d6lTtGZyqIu4/WokHDIftWEYxSc1neZPjwzPqU37nLOxjdGJKR46GCvCyvLD6PgU0Vgyr4J6UyTM44cniCdKl6qTi96Rcbpa6maNM8gHM7nUJfZTF7PYiU+kKWSmb8Mwis/9B2KMJVNz+qd9juZTl4+fuscvHZpXjdobK2hade/weEEDycDz1dZWVZQ8lzoa90zQxfJRA0TCtWb6Ngyj+PhC9+wN8wvqSGOIjWtWldUGHX50dj5Ss3z8euFB8lOrKo+PTBQ0kAy8VL3O5rqZCPNSER1NUFkhtDYUUVA3hhgcS5bM9WOC2jCOU3bsHWJzJEzbAv/gbdvYyp8fGQ7MDkrz0RON01BbxQl5DDpat7qO+prKQAnq0Ykp4slUQQPJfLpaSr8vdTSWYE1DbUFKpc5GpCnEdFoZGiuNVm2C2jCOQyZTae7YvzD/tM85G9s4MjnNvX2jBVxZ/tjVH6c70pBXv21FhbApEg6UoPZ9xoWqSpZJZ3NdIEzfxTR7g2f6htJVJzNBbRjHIXf3HiYxlV6UoPZN5DvLIE1LVfMe8e3T3d7Arv7g+KgLub1lNl3N9cQSqZLWvY6OJooaSAZHA9dKVfTEBLVhHIds33sIETj7pIUL6pZVNZzSES6LgLJDY5OMjE8VRFBv7ghzaCxZMjNoNsUoduLjz1FKrToaL52gLlWKlglqwzgO2bF3iKed0ERTffWi+p2zsY0794+QTAVzYwof3zS9OY+BZD7dAYv87h0Zp6mumsbQ4n7LpdBV4qInialpDo9PFd303dZQQ4WUrjqZCWrDOM6YmJzmr48dXpTZ2+ecja0kU2n++tjhAqwsf8xEfBdIowbYPRAMP3XfyERRzN5w1LxeqsjvQZea1V5kjbqqsoK2htKlaJmgNozjjLseHWFyenH+aZ8zN7RQIcHfn7onGqdlVQ1tDfmv1NUerqWprjowm3MUI4fap6mumnBtVclM377puaPIgho883c0bhq1YRhFYMe+Q1RVCM9c37Lovo2hap6+rokdAfdT74rG2dSe34hvHxGhOxKMmt+q6jTq4ghqEaGzpb5k+1KXoiqZT6Sxlv5RE9SGYRSB7XuH2NLZRENt1ZL6b9vYxt29hxmfTOV5ZflBVdkdHSuIf9qnOxJmV3+85Ht0D8aTJFPpouRQ+5QyRcs3PRfbR+3NGWIgbqZvwzAKTDwxxb19o5yzsW3JY5yzsZWpaeXO/SN5XFn+ODCaYCyZKoh/2mdzR5hYIlXSspKQkZpVJNO3P1ffyERJHlKisQQ1VRU01RU+cC6bSGOI4SOTJQmkNEFtGMcRd+z3dsBaSH3v2di6vpnqSgmsn7qnv3AR3z5HI79La/72i50UK5jMn2tiappDY5NFm9MnGkvQ0Rgq6OYjs+Fr8QMleDgzQW0YxxE79g5RU1XBGU9qXvIY9TVVnNa1OrB+6l1OeHa3Hw+C2t+HurgaNZRmu8toLFESszcc9YsPlCCgzAS1YRxHbN87xBknriZUXbmscbZtbOO+x0eJBWy7R/A06khj7aJzxBeDF1FeW/LI776RCdaEa5f9ey4GP3CtFClaA7Fk0VOzfI5WJztONGoRuVhEdonIHhH51xzHLxeRB0XkXhG5WUSelHHsNyJyWER+ndXnajfm/SJylYhUZxw7T0TuFpEHRORW1xYSkT+LyD2u/eMZ558kIn8Skd0i8iMRqck49hq3tgdE5L8z2j/r2h4Skf8Qj7Cb138dEpEvufP/l4jc59r/ICJPybqeE0VkTETen9VeKSJ/zbx+EfmuiDySMc9prv0UEdkhIskc47zb3asHROQ9Wcfe6e7lAyLy2dl+R6O8ODw+yYMHY8vyT/ts29BKWuHP+4bzsLL80jNQmNKh2WzuKH3kd+/IOF1FDCQDZgLXih1Qpqr0O9N3KShlGdGiC2oRqQS+CrwQeApwabaQAv4KbFXVLcA1QKawuAJ4fY6hrwZOAZ4O1AFvdvOtBr4GvExVnwq82p2fBM5X1VOB04CLReRsd+wzwBdVdRMwArzJjbUJ+BDwLDfWe1z7OcCzgC3A04BnAs9V1biqnua/gEeBn7k5/ltVn+7aPwt8Iet6vghcn+M63w08lKP9Axlz3e3ahoF3AZ/LPFFEngb8E3AmcCrwEndtiMjzgJcDW9w1HtPXKF927htGlWX5p31OP3E1tVUVgfNTT6ddxHcRBHV3JExPdIx0ibY+BCeoi5Sa5bOqtorWVTVFr042lkwxPjldMtN3c301NZUVJSkjWgqN+kxgj6ruU9VJ4Id4gmEGVb1FVf1/BTuBzoxjNwNPeIxV1evUAfw5o8/fAT9T1cfceQPuXVXVrwFY7V4qXpTC+XgPCADfA17hPv8T8FVVHckcC1AgBNQAtW6saOb6nCBsB253fWMZh1e5MfxzXwHsAx7IGqMTeDHw7ezrz4WqDqjqHUC2ffLJwE5VHVfVFHAr8Ep37K3Ap1U1mXWNxyWj41Pc9WjwtMalsGPvIeqqK9nSuXrZY4WqK9m6vpkdAdug47HhcZKpdF73oJ6NzZEwE1PTPH64NDnFqek0Bw4nipqa5eOlaBX3uo+mZpVGoxYR2htrj5tgsnVAb8b3Ptc2G28it2aZE2fyfj3wG9fUDTSLyO9F5C4ReUPGuZUicjcwANykqn8CWoHDToBlr68b6BaRP4rIThG5GEBVdwC3AAfd6wZVzdZ6LwV+pBk5DSLydhHZi6dRv8u1rQI+CHycJ/Il4F+AXBsCf8q5Cr4oIvM9ct4PPEdEWkWkHngR0JVxjec60/+tIvLMXAOIyFtE5E4RuXNwcHCe6cqXf/vFfbzqGztKbuLMB9v3DvHMk1qoqcrPf/ttG1p56GCM4SPFj/6dDd9nXBSN2j0MlMpPfXA0wXRai5qa5dNZgn2p/Trb7eHSCGpw1cmOE406V1x9TtuRiLwO2Ipn7l4oXwNuU9Xb3fcq4Bl4mugLgP8jIt0AqjrtTM+dwJnOJDzX+qqATcB5eIL32yKyWkROxtNSO/GE+vki8pysMS4BfnDMoKpfVdWNeIL5w67543hm92Mq/ovIS4ABVb0rx/o+hGf2fybQ4sabFfcQ8RngJrwHmnsA/8GkCmgGzgY+APxYcuRCqOo3VXWrqm5ds2bNXNOVLQ8djPHrew+iCl/6bU+pl7MsBuNJdg+M5cXs7bPN+bqDtO2l/0B1cntDwefa5ObYVaKHuKPbWxZfUHc113Pg8ATTRTT7z5QPbSqloK49bkzffRzV3sATbgeyTxKRC4F/w/MtL8jWICIfBdYAl2fN9xtVPaKqh4Db8PyyM6jqYeD3wMXAIWC1iPhlmzLX1wf8UlWnVPURYBee4H4lnil5zAnY6/EEnb+uU4GqWYQseOZ/37x+FvBZEdmP5wP/3yLyDjwf+Mtc+w/xHgb+y63/oDPlJ4Hv4LkX5kRVr1TVM1T1OXi+7N0Z1/gzN96f8bT35UcflSFfvKmHcG0Vl217Etfd188DB0ZLvaQl45uot23In6De0tnEqprKQG172RON09VSx6olVl1bDOFQNetW15XM2uJHXZdCo+5qqWNqWouqXfqm7/ZwaXzU4KqTHSem7zuATS6yugZP07w28wQROR34TzwhvSAfqYi8GU9jvlRVM03Dv8Qz5VY5M+9ZwEMissYFmiEidcCFwMPONH0L8CrX/zI3BsAvgOe5Pm14ZuJ9wGPAc90c1cBzOTbg61KytGk/eMvxYpygVNVzVXW9qq7HM3X/u6p+RVU/pKqdrv0S4Heq+jo31lr3LngC//4F3K92934i8DcZ6/sFno8eZ3mowXt4Oa64t+8wNz4Y5c3nbuDyizbTGKriizeVr1a9Y+8hwqEqnnpCY97GrK6s4MyTWtgRoICynmi8KGZvn+5IQ8lM333D41QIrF1dfA1zJpe6iJHf0ViCcG1VUR7CZiPSGGIsmWIsWdzyuUUX1M73+w7gBjxh9mNVfUBEPiEiL3OnXQE0AD9x6UYzglxEbgd+AlwgIn0i8gJ36BtABNjh+nzEzfcQnnn3Xrwgs2+r6v3AWuAWEbkX7+HhJlX1U54+CFwuInvwfNZXuvYbgCEReRBPmH9AVYfwAs/2AvfhmZHvUdVfZVz2a8gS1MA7XPrT3XgWgMsWey8zuFpE7nPztwH/192rDhHpc+N/2N0v/y/1T911/Ap4ux8gB1wFbBCR+/E098sy/erHC1+4qYfV9dX847PX01RXzVues4HfPjTA3b3B3t5xNnbsHeKsk1qpqszvf/ltG1vZO3ikJH67bCZTafYNHilKapZPd0eYfYNHSE3nChspLL0jE6xtqqM6z7/pQvDN7cXcnCMaSxApodkbMquTFfffe0keTVT1OuC6rLaPZHy+cI6+587SPuu1qOoVZPm5VfVe4PRZzt9HDvOxE1iXc6xpHVWdBv55jvk35Gh792znZ5zzsVnaf49nqve/nz/Lef1kRMxnHZvtPk4Cr5tvbSuZux4d5ve7BvngxacQDnnp+G981klc+YdH+MJNPXz/H+f1LASKxw9PsH9onNdvW5/3sf2c7B17h3jF6XPFhBaeRw4dIZXWgpYOzWZzJMzkdJr9Q+NF8Ytn0js8XpKIb4ATVocQKb5GXarULJ+IC2TrjyXYsKZ4v7dVJjOMLD5/Yw9tDTVcds5MnR0aaqv4X8/dyG09g9yxv7zStXzTdD4DyXyevLaRprrqQJi//aCuTQUsHZpNKUuJliKH2qe2qpJIOFTUyO9oLDkjKEuFr9EX209tgtowMtixd4jte4d463knU19zrJHmDdvW09ZQy+dv3FWi1S2NHXuHaFlVUxDfbWWFcNZJLWzfV/owht3ROJUVwoY1q4o258ntDYgUP0UrMTVNNJYsSSCZT1dLHX1FyqVWVQbiiZKVD/UpVXUyE9SG4VBVvnDTLiKNtfz9WSc+4XhdTSVvf95Gdu4bZvue0gumhaCq7Nh7iLM3tFBRUZgdh87Z2Erv8ETJ9ij22dUfZ31rfVHrXoeqK1nfuordA8UV1H6RlWLumpVNV3PxcqmHj0wyNa10lNj03VBbxaqayqLX+zZBbRiO23Yf4o79I7zj/E2z/rG/9MwTWdsU4nM37irJfryL5dGhcQ6MJmZyngvBOSc7P3WJ86l7ovGi+qd9ShH5PZOaVSLTN3hFT/pjCSZThQ+kK3VVskxKUfTEBLVh4LTpG3exbnUdr93aNet5oepK3nH+yfzlscP8vif4Fdm2F9A/7bOpvYG2hpqS+qknJqd5dHi8qBHfPpsjYfYPjZOYmi7anL71oqSm7+Y6VOFAEUqoRt3WkqU2fYMJasMoGb99aIB7+kZ51wUnz1ti89XP6KKzuY4v3NgTeK16x74hIo21bGgrnN9WRDh7Qyvb9x4q2f3YMzCGKiUR1JsiYabTyr7BI0Wbs3dknJrKipIW/ziaolV483d0tPRVyXwijbUzDw7FwgS1cdyTTitfuKmH9a31/M0ZObPZjqGmqoJ3X7CJ+x4f5cYHo/OeXyp8//S2Da3kqAKbV87Z2EY0lmTfoeIJq0z8qOuSaNQdxY/87hueYF1zXcHiDhbC0e0ui6BRO9P3mobS+qjB16iTRX0oNUFtHPf85oF+HjoY490Xblpw8YhXnr6ODW2r+OJNPSXd5nAudg+McWhsMi/7T8+Hb1ovlfm7JxqnprKC9a3FNwWvb11FdaUUVVD3jpQuh9pnbVMdVRVSHI06nqB1VU3eNpRZDpHGEJOpNIfHszclLBylv2rDKCHTTps+ub2Bl5268IIdVZUVvPvCTTzcH+d/7jtYwBUuHV9obiugf9rnSa31rG0KlUxQ74rG2djekPfKawuhpqqCDW0NxRXUw6XLofaprBBOWF1XlGj/gVjpU7N8ZlK0imj+NkFtHNf86p4D7BkY470XdlO5SDPiS7ecQGOeUTkAABQKSURBVHekgS/+tqckJSTnY/veQ3Q21xXlD7qIsG1jKzv2DZXEwtDTH2dzpLiVwTLp7ggXbRetsWSKkfGpkgaS+XS11M1EoBeS/lii5KlZPn51tGKmaJmgNo5bUtNpvvTbHk7pCPPCp3Usun9FhfDeC7vZN3iEX979hA3gSko6rezcN1zQaO9sztnYxvCRSXqKnFMcS0xxYDTBphL4p3262xvoHZ7gSBE2a+ib2d6ytKZv8KLO+4ph+o4lA5GaBRka9ahp1IZRcH72l8fZPzTO+y7avOSgnBc8tYOnntDIl2/ezVSAtOoHD8YYnZgqin/axzexb99TXPP37qi3dXsxd83KptsFlO0eGJvnzOXjB28FQ6Ou59DYJOOThXtASU2nOTSWDIzpu31GozZBbRgFZTKV5ss37+bUziYufHL7ksepqBAuf343jw2P89O7+vK4wuVRTP+0z7rVdTyptX4md7tY+L7hUhQ78dlcxJrfvk+41MFkmWsopPl7cCyJKnQERFDXVlXSXF9tPmrDKDQ/urOXxw9PcPlFm5edunT+Ke2c1rWa/7h5N8lU8YpezMX2vYfYsGZV0c2F52xs5U+PDDFdRD/1rv449TWVrFtdwnKaLfWEqivoKUKFst6RceprKmlZVVPwueajswj7Uh+tShYMHzV45u/+UfNRG0bBSExN89Xf7WHrk5p5zqblm4ZFhPdd1M2B0QQ/uqM3DytcHlPTaf78SHH90z5nb2glnkjxwIHRos3ZE42zKRIuaU5xZYWwqb04AWW9wxN0NdcXPDd+Ifh+8sIKak9zDYqPGry1DJhGbRiF47//9Bj9sQSXX9Sdtz92zz65jTPXt/CV3+0painJXNz3+ChHJqfZtqF4/mmfGT91Ec3fPdE43UXeCzoXmyLFSdHqGxkPRCAZeAVIQtUV9BbQ9D0Q88uHBkmjrjUftWEUivHJFF/7/R7O2dia10ArX6seiCf5r52P5m3cpeD7p8/e0FL0udvDITa1NxRNUA+NJTk0NllS/7TP5kiYaCzJaAELYagqfSMTMybnUiMidBY48rs/lqCyQmhbFSRBHWIwniyai8cEtXFc8f0dj3JobJL3XdSd97HP2tDKs09u4+u/31uUNJ3Z2LF3iFM6wrSWqNziORtbuXP/cFF2VepxEd+lKB2ajR/5Xcj0tMPjU4wlUyUvdpJJV3NdQcuIRmNJ2sO1JXVtZNPeGCKtcGisOH5qE9TGcUM8McU3bt3LeZvX8IwnFUbbvPyiboaOTPLd7fsLMv58JFPT3LF/uKjR3tls29jK+OQ09/YdLvhcQYj49vEjvwu55aVfrrMrABHfPl0thd2XOhqgqmQ+fgR6sczfJqiN44bv/HE/h8enuPz5+demfc44sZnzT2nnm7ftI5YoXi1gn78+dphkKl3U/OlszjqpFZHi+Kl3ReM0hqpKuouUz9qmEOHaqoL6qX3NNSimb/BStOKJVMFM/gOxJJEA/L6ZFLs6mQlq47hgdHyKb92+j+c/JcKWztUFnevy53czOjHFVX94pKDz5GL73iEqBM48qfj+aZ/mVTU8ZW1jUep+9/TH2dwRDkQEtIiwKdJQHI06IMFkcLTwSqG06v5YIhDbW2YSMY3aMPLPt27fRzyRKqg27fO0dU1c/NQOrrz9EQ6PTxZ8vkx27h3i6euaaKqrLuq82Wzb0Mpdj40UNAJeVb2I7wD4p302d4TpicYLtgVi7/A4q+urCYdK+/tmMrMvdQFStBJT04xOTAUqNQugdVUNFWKC2jDyxtBYku/88RFevGUtT17bWJQ53/v8bsYmU3zztn1FmQ+8iPa/9o5wdgn90z7nnNzKZCrNXx4dKdgc0ViSWCIVCP+0T3ckzMj4FIfGCvOA1jsyEYjSoZkUUqMecKblILg2MqmqrGBNuHgpWiaojRXPf962j4mpad574aaizbm5I8xLt5zAd/64v2iRoXfuH2FqWkvqn/Z55voWKiuEHfsKZ/72i4sESqMucCnRIOVQ+zTVVxMOVRWkjGi/E4RBM32DZ/42H7Vh5IGBeILv79jPK05bx8ntxf2D/u4LN5FMTfON3+8tynw79g1R9f+3d+9BVtb3Hcffn91lue3CAsvFwMICAtZxFBSVi1EUbUnqSDIJkxq1tGYGm8ZLMs1YbTo2kz9SUjsGp9rYNDUyhpCk1KSYJqIi6EREQeRmkJs3Vq7CclmW6+63fzy/A4fDLtfn7PMDvq+ZM/uc5zz7e76HZfd7fr/n9/y+JeLq2m5tcr4TqezQjsv7dS3qhLLccp0xJeohRZz53dyc3EMdW48akl51MYa+Y1yVLKdXZQfvUTuXhn+ft55DTcb949uuN50zuGcFXxzRj2cXftQmv9AL1m9neE0VncrLin6uUzF6UA+WbdhJQ5HuKV+9ZQ/VFe2jWPM6p7qinO6dy4vSo97WcICDh5ujKMZRqF+3jkVZnexIoq6ML1H36epD386dtY079/HzNz9m0lX9qK3unEkMD4wfQlOz8eS8dUU9z+79h1hRtzOT9b1bM2ZwNYebjUUf7ihK+2u27GFYn+yXDs0niaG9K4qy5veRqlkRLXaSU9M9WZ0s7Ul0W/ccoH1ZCV06xvHhM1/vyg7UNx5qk0I8nqjdeeuJeeswjHtvujizGPr36MSkkTXMfOvjoi6zuOiDHTQbjI7g+nTOVQO6UV5awsIiDH83NxtrtzRENeydM6x3JWu3NKSetI4udhJhou7Wkf2HmtmW8nyMzbuSW7NiuP2uUG44fmsbXKf2RO3OSxt2NPKrRRu4/Zr+mS8Ocd9NFyPEE68Ur1e9YP12ystKGNG/uPeIn46O5aUM719VlOvUdfX72Heo6cjkrZgM7VNJw4HDbNyV7rDo0cVO4hv6PnqLVrrD31t2749y2BuOFglpi+FvT9TuvPT43LWUlohv3JhdbzrnM1Ud+eq1/fnvt+v4aPveopxjwfrtjBzQjQ7tSovS/pkaM7gHKzfuSn3VqtzQ8pAYE3Vu5nfKE8rq6hvpVdk+up8xHE3UaY8abd1zIKqqWflyM9HbYua3J+oLjKQJklZLWifpoRZeby/pl+H1NyXVtn2UZ+f9bQ08t6SOO0cNiGa26N+OG0xZiXh87trU267fe5BVm3ZHdX06Z8zgaszgzQ/S7VWvOXJrVlzXqAGGhrsL0r5OvWHHvqiKceTL9fLTvEXLzJKh70h+hwvlevreo3apklQKPAl8DrgUuF3SpQWHfQ2oN7OLgR8CP2jbKM/etJfX0r6slK+PG5x1KEf06tKByWNq+c07n7Bua0OqbS8M9ypnWYijNVfUdKVDu5LUh79Xb95D36qOUa3QldO1Uzv6dOmQeo96Q31jVMU48nUqL6O6ojzVW7T2HDjMvkNN0XzYLlTVqR3lZSWeqF3qrgHWmdn7ZnYQ+AUwseCYicD0sD0LGK8YZ3K0YvXmPTy/fCN/NbaW6ozKPLbmnusH0aFdKdNeXpNquwvWb6dTeWnR1zA/E+3LSrm6tnvq634nS4fG15vOGdqnMtVyl4ebmtm0a3/m8y1OpG+3dKtobQ0JMNahb0n07tI2t2jFN+fdFVNfYEPe8zrg2taOMbPDknYBPYBP0w5mZ+NBJj31Rqpt1jcepKK8jHuuH5Rqu2noUdGeu8cO5Il561i9+dXU2q2r38e1g7rTrjTOz92jB/fgX15YzS2Ppfee129r4IZhPVNrL23Delfw9Oufpvaem5qNpmaLblWyfDXdOvLiH7ek9p73hXXiY+1RQzL83RbXqD1RX1ha6hkX3kNyKscgaQowBaB///5nFExJSVJtKG0TLruIqk7xLIKRb8oNg9iyez97D6a3CMjQ3pXcNXpAau2l7Ysj+rJm8x4ONjWn1uYlF3XhS1f2S629tH35qho27dpPc4q3aF1RU8WNw3ql1l7a7hw1ADOw4/9cnLGxg6sZXhPfSFHOuGE9aTxY/PuoVawqLy4+kkYD3zWzPwvPHwYws3/OO2ZOOOYNSWXAZqCnneA/ysiRI23x4sXFDd45584zkt42s5EnOy7OsTJXLIuAIZIGSioH/gKYXXDMbGBy2P4y8MqJkrRzzrni8qHvC0i45nwvMAcoBZ42s3clfQ9YbGazgf8CnpW0DthBksydc85lxBP1BcbMfgf8rmDfI3nb+4FJbR2Xc865lvnQt3POORcxT9TOOedcxDxRO+eccxHzRO2cc85FzBO1c845FzFf8MSdNUnbgI/OoolqirBEaYpijw/ijzH2+MBjTEPs8UFcMQ4ws5OuheuJ2mVO0uJTWZ0nK7HHB/HHGHt84DGmIfb44NyIsZAPfTvnnHMR80TtnHPORcwTtYvBj7MO4CRijw/ijzH2+MBjTEPs8cG5EeMx/Bq1c845FzHvUTvnnHMR80TtMiNpgqTVktZJeijreApJqpE0T9IqSe9KeiDrmFoiqVTSO5J+m3UsLZFUJWmWpPfCv+XorGPKJ+lb4ee7UtJMSR0iiOlpSVslrczb113SS5LWhq/dIozx0fBzXi7p15KqYosx77VvSzJJ1VnEdjo8UbtMSCoFngQ+B1wK3C7p0myjOs5h4O/M7E+AUcA3IowR4AFgVdZBnMDjwAtmdglwBRHFKqkvcD8w0swuIyn/GkNp12eACQX7HgLmmtkQYG54nqVnOD7Gl4DLzOxyYA3wcFsHVeAZjo8RSTXALcDHbR3QmfBE7bJyDbDOzN43s4PAL4CJGcd0DDPbZGZLwvYekgTTN9uojiWpH/DnwE+yjqUlkroA15PUOcfMDprZzmyjOk4Z0FFSGdAJ2JhxPJjZayT14PNNBKaH7enAF9o0qAItxWhmL5rZ4fB0IdCvzQM7Np6W/h0Bfgg8CJwTk7Q8Ubus9AU25D2vI7IkmE9SLTACeDPbSI4zjeQPTnPWgbRiELAN+GkYnv+JpM5ZB5VjZp8A/0rSs9oE7DKzF7ONqlW9zWwTJB8igV4Zx3MydwO/zzqIQpJuAz4xs2VZx3KqPFG7rKiFfVF+upVUAfwP8E0z2511PDmSbgW2mtnbWcdyAmXAlcCPzGwEsJfsh2yPCNd5JwIDgc8AnSXdmW1U5z5J3yG5dDQj61jySeoEfAd4JOtYTocnapeVOqAm73k/IhhyLCSpHUmSnmFmz2UdT4GxwG2SPiS5dHCTpJ9lG9Jx6oA6M8uNRMwiSdyxuBn4wMy2mdkh4DlgTMYxtWaLpIsAwtetGcfTIkmTgVuBOyy++38Hk3woWxZ+b/oBSyT1yTSqk/BE7bKyCBgiaaCkcpIJPLMzjukYkkRybXWVmT2WdTyFzOxhM+tnZrUk/36vmFlUvUEz2wxskDQs7BoP/DHDkAp9DIyS1Cn8vMcT0WS3ArOByWF7MvC/GcbSIkkTgL8HbjOzxqzjKWRmK8ysl5nVht+bOuDK8P80Wp6oXSbChJN7gTkkfxh/ZWbvZhvVccYCd5H0VJeGx+ezDuocdB8wQ9JyYDjw/YzjOSL09GcBS4AVJH8TM1+5StJM4A1gmKQ6SV8DpgK3SFpLMmN5aoQxPgFUAi+F35enIozxnOMrkznnnHMR8x61c845FzFP1M4551zEPFE755xzEfNE7ZxzzkXME7VzzjkXMU/UzrnzjqSGVvYvCF9rJX015XP+Q0vncu5s+e1ZzrmikFRqZk0ZnbvBzCpO8Po44NtmdutptHnC93Oyczp3prxH7Zw7LaE3+p6k6aHu8KywhjKSPpT0iKQ/AJMkDZe0MK8+cbdw3HxJ0yQtCHWgrwn7u0v6TTh+oaTLw/4b8hadeUdSpaQKSXMlLZG0QtJJq6/l9bSnAp8N7X1LSU3vRyUtCue+Jxw/TklN8p+TLIhCiO9tJTWsp4R9U0kqcC2VNCP/XEo8Gt7nCklfyWt7vo7W6p4RVkdz7lhm5g9/+MMfp/wAakkKqIwNz58m6Z0CfAg8mHfscuCGsP09YFrYng/8Z9i+HlgZtv8N+KewfROwNGw/n3e+CpJiH2VAl7CvGljH0VHChlZibwhfxwG/zds/BfjHsN0eWEyyJvQ4kkIiA/OO7R6+dgRWAj1aOmfeub5EUqe5FOhNsmzpRaHtXSTrTZeQrKB1XdY/X3/E9/AetXPuTGwws9fD9s+A6/Je+yWApK5AlZm9GvZPJ0nKOTPhSM3gLpKqQjvPhv2vAD1CO68Dj0m6P7R5mKQC2/fD0qQvk5RJ7X2G7+dPgb+UtJSklGkPYEh47S0z+yDv2PslLSOpt1yTd1xrrgNmmlmTmW0BXgWuzmu7zsyagaUkH4KcO0ZZ1gE4585JhZNb8p/vPYs2Wix/amZTJf0f8HlgoaSbgVFAT+AqMzsUqiF1OMVzFxJwn5nNOWZnci17b8Hzm4HRZtYoaf4pnPNEw9kH8rab8L/JrgXeo3bOnYn+kkaH7duBPxQeYGa7gHpJnw277iLpTebkrtVeB+wKx78G3BH2jwM+NbPdkgZbUvnoByTD0pcAXUnqcR+SdCMw4DTi30NSPCJnDvB1JWVNkTRUUucWvq8rUB+S9CUkHxZyDuW+v8BrwFfCdfCeJKMKb51GrO4C55/enHNnYhUwWdJ/AGuBH7Vy3GTgqTDZ7H3gr/Neqw+3MHUB7g77vgv8NAxnN3K0rOM3QzJuIimT+XuSRPu8pMUkw8bvnUb8y4HDYQj7GeBxkmHnJWFC1zbgCy183wvA34T4VpMMf+f8GFguaYmZ3ZG3/9fAaGAZyajBg2a2OSR6507Kb89yzp0WSbUkE7EuO4s25pNMQFucUljOnbd86Ns555yLmPeonXPOuYh5j9o555yLmCdq55xzLmKeqJ1zzrmIeaJ2zjnnIuaJ2jnnnIuYJ2rnnHMuYv8P4tijnG0/OXAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from io import StringIO\n",
    "import sys\n",
    "\n",
    "x = np.random.randn(100)\n",
    "\n",
    "old_stdout = sys.stdout\n",
    "mystdout = sys.stdout = StringIO()\n",
    "\n",
    "with pm.Model() as model:\n",
    "    mu = pm.Normal('mu', mu=0, sigma=1)\n",
    "    sd = pm.Normal('sd', mu=0, sigma=1)\n",
    "\n",
    "    mu_print = tt.printing.Print('mu')(mu)\n",
    "    sd_print = tt.printing.Print('sd')(sd)\n",
    "\n",
    "    obs = pm.Normal('obs', mu=mu_print, sigma=sd_print, observed=x)\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(5, step, tune=0, chains=1, progressbar=False) # Make sure not to draw too many samples\n",
    "\n",
    "sys.stdout = old_stdout\n",
    "\n",
    "output = mystdout.getvalue().split('\\n')\n",
    "mulines = [s for s in output if 'mu' in s]\n",
    "\n",
    "muvals = [line.split()[-1] for line in mulines]\n",
    "plt.plot(np.arange(0, len(muvals)), muvals)\n",
    "plt.xlabel('proposal iteration')\n",
    "plt.ylabel('mu value');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['mu']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that for each iteration, 3 values were printed and recorded. The printed values are the original value (last sample), the proposed value and the accepted value. Plus the starting value in the very beginning, we recorded in total `1+3*5=16` value above."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
